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Abstract. 

We study the scaling of the Renyi entanglement entropy of two disjoint blocks 
of critical lattice models described by conformal field theories with central charge 
c = 1. We provide the analytic conformal field theory result for the second order 
Renyi entropy for a free boson compactified on an orbifold describing the scaling 
limit of the Ashkin- Teller (AT) model on the self-dual line. We have checked this 
prediction in cluster Monte Carlo simulations of the classical two dimensional AT 
model. We have also performed extensive numerical simulations of the anisotropic 
Heisenberg quantum spin-chain with tree-tensor network techniques that allowed 
to obtain the reduced density matrices of disjoint blocks of the spin-chain and to 
check the correctness of the predictions for Renyi and entanglement entropies from 
conformal field theory. In order to match these predictions, we have extrapolated 
the numerical results by properly taking into account the corrections induced by 
the finite length of the blocks to the leading scaling behavior. 



Entanglement of 2 disjoint intervals in c — 1 theories 



2 



1. Introduction 

Let us imagine to divide the Hilbert space H of a given quantum system into two parts 
Ha and Hb such that V. = Ha ® Hb- When the system is in a pure state |*), the 
bipartite entanglement between A and its complement B, can be measured in terms 
of the Renyi entropies [JJ 

S^ = -^—]ogTr^, (1) 
1 — n 

where pa — Tr^ p is the reduced density matrix of the subsystem A, and p = |*)( V E'| 

(n) 

is the density matrix of the whole system. The knowledge of S A as a function of n 
identifies univocally the full spectrum of non-zero eigenvalues of pa [2], and provides 
complementary information about the entanglement to the one obtained from the von 
Neumann entanglement entropy ■ Furthermore, the scaling of with the size 
of A in the ground-state of a one-dimensional system is more suited than to 
understand if a faithful representation of the state in term of a matrix product state 
can be or cannot be obtained with polynomial resources in the length of the chain 

mm- 

For a one-dimensional critical system whose scaling limit is described by a 
conformal field theory (CFT), in the case when A is an interval of length I embedded 
in an infinite system, the asymptotic large t behavior of the quantities determining 
the Renyi entropies is El [7j |8] 

/ ,x c(n-l/n)/6 , x x £ 

Trp'X^cJ-) , ^sW-Efl + AUgi+c;, (2) 

\a J o \ n J a 

where c is the central charge of the underlying CFT and a the inverse of an ultraviolet 
cutoff (e.g. the lattice spacing). The prefactors c„ (and so the additive constants c' n ) 
are non universal constants (that however satisfy universal relations [5]). 

The central charge is an ubiquitous and fundamental feature of a conformal field 
theory [TO], but it does not always identify the universality class of the theory. A 
relevant class of relativistic massless quantum field theories are the c = 1 models, 
which describe many physical systems of experimental and theoretical interest. The 
one-dimensional Bose gas with repulsive interaction, the (anisotropic) Heisenberg 
spin chains, the Ashkin- Teller model and many others are all described (in their 
gapless phases) by c = 1 theories. These are all free-bosonic field theories where 
the boson field satisfies different periodicity constraints, i.e. it is compactificd on 
a specific target space. The two most notable examples are the compactification 
on a circle (corresponding to the Luttinger liquid field theory) and on a Z-i orbifold 
(corresponding to the Ashkin- Teller model [TTJ Q21 [IB] ) . The critical exponents depend 
in a continuous way on the compactification radius of the bosonic field. A survey 
of the CFTs compactified on a circle or on a Z^ orbifold is given in Fig. [T] in a 
standard representation [121 113] . The horizontal axis is the compactification radius 
on the circle r c i rc i c , while the vertical axis represents the value of the Z^ orbifold 
compactification radius r or b. The two axes cross in a single point, meaning that the 
theories at r C i rc i c = y/2 and at r or b = 1/ \/2 are the same. (The graph is not a cartesian 
plot, i.e. it has no meaning to have one r C i rc i e and one r or b at the same time.) For 
some values of r c i rc i c and r or b, we report statistical mechanical models and/or field 
theories to which they correspond. In the following we will consider the Ashkin- Teller 
model that on the self-dual line is described by r or b £ [-\/2/3, a/2] and the XXZ spin 
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Figure 1. Survey of c = 1 theories corresponding to a free boson compactified 
on a circle (horizontal axis) and on an orbifold (vertical axis) as reported e.g. 
in Refs. 1 1 21 . For some values of r c ; rc i c and r oz \,, the corresponding statistical 
mechanical models are reported. The XXZ spin chain in zero magnetic field lies 
on the horizontal axis in the interval r c ; rc i e S [0, l/\/2]. The self-dual line of the 
Ashkin- Teller model lies on the vertical axis in the interval r or t, £ [*/2/3, v2]. 



chain in zero magnetic field that is described by r c ; rc i c £ [0, l/v2]- We mention that 
different compactifications have been studied [TJ] , but they correspond to more exotic 
statistical mechanical models and will not be considered here. 

According to Eq. ([2]), the central charge of the CFT can be extracted from the 
scaling of both the Renyi and von Neumann entropies. In the last years, this idea has 
overcome the previously available techniques of determining c, e.g. by measuring the 
finite size corrections to the ground state energy of a spin chain [15] . However, the 
dependence of the scaling of the entropies of a single block only on the central charge 
prevents to extract from them other important parameter of the model such as the 
compactification radius. It has been shown that instead the entanglement entropies of 
disjoint intervals are sensitive to the full operator content of the CFT and in particular 
they depend on the compactification radius and on the symmetries of the target space. 
Thus they encode complementary information about the underlying conformal field 
theory of a given critical quantum/statistical system to the knowledge of the central 
charge present in the scaling of the single block entropies. (Oppositely in 2D systems 
with conformal invariant wave-function, the entanglement entropy of a single region 
depends on the compactification radius US].) 

This observation boosted an intense theoretical activity aimed at determining 
Renyi entropies of disjoint intervals both analytically and numerically [TTJ HH [TU [201 
iniE2[23l[^[25l[26l|2Il[Ml[Ml[3Q]. A part of this paper is dedicated to consolidate 
some of the results already provided in other works where they either have been studied 
only on very small chains, with the impossibility of properly taking into account the 
severe finite size corrections [T7] or have been tested in the specific cases of spin chains 
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equivalent to free fcrmionic models 24, 26] . An important point to recall when dealing 
with more than one interval is that the Renyi entropies in Eq. ([!]) measure only the 
entanglement of the disjoint intervals with the rest of the system. They do not measure 
the entanglement of one interval with respect to the other, that instead requires the 
definition of more complicated quantities because A% U A 2 is in a mixed state (see e.g. 
Refs. [31 for a discussion of this and examples). Furthermore, it must be mentioned 
that some results about the entanglement of two disjoint intervals are at the basis of 
a recent proposal to "measure" the entanglement entropy [32] . 



1.1. Summary of some CFT results for the entanglement of two disjoint intervals 

We consider the case of two disjoint intervals A = A\ U A 2 = [tti,«i] U [1*2 > 1*2] • By 
global conformal invariance, in the thermodynamic limit, Tr//| can be written as 

\(n-l/n) 



ui - u 2 v 1 



V 2 



F n {x) 



(3) 



(4) 



\ui - v x \\u 2 - vi\\u\ - v 2 \\u 2 - v\ 

where x is the four-point ratio (for real uj and Vj, x is real) 

(ui - Vi)(u 2 - v 2 ) 

x —— 

(ui - u 2 )(vi - v 2 ) ' 

The function F n (x) is a universal function (after being normalized such that ^(O) = 1) 
that encodes all the information about the operator spectrum of the CFT and in 
particular about the compactification radius. c„ is the same non-universal constant 
appearing in Eq. 

Furukawa, Pasquier, and Shiraishi [T7] calculated F 2 (x) for a free boson 
compactified on a circle of radius r C j rc i c 

9 3 ( V r)e 3 (T/ v ) 
F2{x) = WjW ' 

where Q v are Jacobi theta functions and the (pure-imaginary) r is given by 

4 



(5) 



i(r) 



t(x) 



2 F 1 (l/2,l/2;l;l-a 
2^(1/2,1/2;!; x) 



(6) 



77 is a universal critical exponent related to the compactification radius n = 2r^ irc | c . [|] 
This has been extended to general integers n > 2 in Ref. [TH] 

9( 01^)6(0^/7?) 

[e(o|r)]2 ' 

where T is an (n — 1) X (n — 1) matrix with elements [19 

n-l 



F n (x) = 



(7) 



r — 



2i 



fc = i 



and 



sin f 71 "^ Pk In cos 

2 Fi(y, l-y;l;l-x) 
2F 1 (y,l-y;l;x) 



27T — (r — s) 
n 



(8) 



(9) 



X Because of the symmetry r) — > 1/r] or r c i rc i c — > l/2r c i rc i e for any conformal property one could 
also define r\ = l/2r^ irclo as sometimes done in the literature. However, corrections to scaling are 
not symmetric in rj — > 1/r) and this is often source of confusion. A lot of care should be used when 
referring to one or another notation. 
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r\ is the same as above, while G is the Riemann-Siegel theta function 
0(O|r) = exp^Trm'-r-m] . 



(10) 



m e z° 



The analytic continuation of Eq. |7| to real n for general values of n and x (to obtain 
the von Neumann entanglement entropy) is still an open problem, but results for 
i«l and 77 <C 1 are analytically known [HI [30] . 

The function F n (x) is known exactly for arbitrary integral n also for the critical 
Ising field theory 30J. However, in the following we will need it only at n = 2 (i.e. 
F2(x)) for which it assumes the simple form [53] 

1 1/2 



Fl s (x) 



1/2 



+ x 1 ' i + {{l-x)x) 1 l i + {l-x) 1 l i 



.(11) 



1 

In Ref. [SO], it has been proved that in any CFT the function F n (x) admits the 
small x expansion 



F„ 



1 + 



( x 



S2{n) 



2a 



s 4 (n) 



(12) 



Un 2 / 

where a is the lowest scaling dimension of the theory. The functions Sj(n) are 
calculable from a modification of the short-distance expansion [30] . and in particular 
it has been found 130 



S2[n) 



n-l 

3=1 



1 



2a 



(13) 



where the integer M counts the number of inequivalent correlation functions giving 
the same contribution. This expansion has been tested against the exact results for 
the free compactified boson (Ising model) with a = minfrj, I/77] (a = 1/4) and JV = 2 
CAf = l). 

All the results we reported so far are valid for an infinite system. Numerical 
simulations are instead performed for finite, but large, system sizes. According to 
CFT [S], we obtain the correct result for a chain of finite length L by replacing all 
distances «y with the chord distance L/n sm(nUij/L) (but different finite size forms 
exist for excited states [23]). In particular the single interval entanglement is [S] 

-c(n—X/n)/6 



L 



sin 



TTfl 



(14) 



and for two intervals, in the case the two subsystems A\ and Ai have the same length 
£ and are placed at distance r, the four-point ratio x is 

2 



sin n£/L 
sin7r(£ + r)/L 



(15) 



1.2. Organization of the paper 

In this paper we provide accurate numerical tests for the functions F n (x) in truly 
interacting lattice models described by a CFT with c = 1. In Sec. [2] we derive the CFT 
prediction for the function F^{x) of a free boson compactified on an orbifold describing, 
among the other things, the self-dual line of the AT model when r or b G [y/2/3, \/2]. 
In order to check this result, we needed to develop a classical Monte Carlo algorithm 
in Sec. [3] based on the ideas introduced in Ref. [TS]. This algorithm is used in Sec. [i] 
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to determine F2{x) for several points on the self-dual line. We also consider the XXZ 
spin-chain in zero magnetic field to test the correctness of Eq. ([7]) . In order to extend 
the results of Ref. |17j to longer chains, we have used a tree tensor network algorithm 
that has allowed us to study chains of length up to L = 128 with periodic boundary 
conditions. In this way we have been able to perform a detailed finite size analysis that 
was difficult solely with the data from exact diagonalization reported in Ref. [17] . The 
analysis also shows that only through the knowledge of the unusual corrections to the 
leading scaling behavior [331 [3D] [37[ [3SJ [2D] we are able to perform a quantitative 
test of Eq. ([7|. The tree tensor network algorithm is described in Sec. [5] while the 
numerical results are presented in Sec. [6] The various sections are independent one 
from each other, so that readers interested only in some results should have an easy 
access to them without reading the whole paper. 

2. n = 2 Renyi entanglement entropy for two intervals in the 
Ashkin- Teller model 

In a quantum field theory Tr p\ for integer n is proportional to the partition function 
on an n-sheeted Riemann surface with branch cuts along the subsystem A, i.e. 
Trp^ = Z n (A)/Z™ where Z n (A) is the partition function of the field theory on a 
conifold where n copies of the manifold 1Z = system x R 1 are coupled along branch 
cuts along each connected piece of A at a time-slice t — [81 139] ■ Specializing to CFT, 
for a single interval on the infinite line, this equivalence leads to Eq. ^ 6], whose 
analytic continuation to non-integer n is straightforward. When the subsystem A 
consists of N disjoint intervals (always in an infinite system), the n-sheeted Riemann 
surface lZ n .N has genus (n — 1)(N — 1) and cannot be mapped to the complex plane 
so that the CFT calculations become more complicated. 

However, for two intervals (N = 2), when for a given theory the partition function 
on a generic Riemann surface of genus g with arbitrary period matrix is known, Tr p\ 
can be easily deduced exploiting the results of Refs. [TH1 [3D]. In fact, a by-product 
of the calculation for the free boson [TDJ is that the (n — 1) x (n — 1) period matrix 
is always given by Eq. Although derived for a free boson, the period matrix is 
a pure geometrical object and it is only related to the structure of the world-sheet 
1Z 71 ,2 and so it is the same for any theory. This property has been used in Ref. [30] to 
obtain F n (x) for the Ising universality class for any n, in agreement with previously 
known numerical results [2DJ. When also n = 2, the surface 7vL 2 ,2 is topologically 
equivalent to a torus for which the partition function is known for most of the CFT. 
The torus modular parameter r is related to the four-point ratio by Eq. Thus, 
the function F2(x) is proportional to the torus partition function where r is given by 
Eq. (|6| and with the proportionality constant fixed by requiring F 2 (0) = 1. This way 
of calculating is much easier than the general one for gDJ HD] and indeed it 
has been used to obtain the first results both for the free compactificd boson [T7] and 
for the Ising model [23] . 

For a conformal free bosonic theory with action 

S= - 1 - / dzdz dcj)d<t> , (16) 

the torus partition functions are known exactly both for circle and orbifold 
compactification [4H l42l IT2] . 
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We now recall some well-known facts in order to fix the notations and derive 
the function F 2 (x) for the Ashkin- Teller model. The bosonic field <j> is said to be 
compactified on a circle of radius r c i rc i c when eft — eft + 2nr c i rc \ c . The torus partition 
function (and the one on the n-sheeted Riemann surface) should be derived with this 
constraint. It is a standard CFT exercise to calculate the resulting torus partition 
function gH [12] 

O 3 (t)t)0 3 (t/t,) 

-^circle 7?) = ! / xTn > ( 17 ) 

where T)d(t) is the Dedekind eta function and r/ — 2rj? irclc . Using Eq. ^ and some 
properties of the elliptic functions, Eq. ([5| for F 2 (x) follows [TTj . When specialized 
at 77 = 1/2 (or i] = 2), F 2 (x) has the simple form 

(a:) = > /(l+ a: V2)(l + (l- a ;)i/a)/2, (18) 

that describes the XX spin-chain (that is equivalent to free fermions via the non-local 
Jordan- Wigner transformation) . 

The concept of orbifold emerges naturally in the context of theories whose Hilbert 
space admits some discrete symmetries. Let us assume that G is a discrete symmetry. 
For the free bosonic theory, the simplest example is the one we are interested in, i.e. 
the Z 2 symmetry. It acts on the point of the circle S 1 in the following way 

g:<t>^-<ft- (19) 
For the partition function of a theory on the torus, we introduce the notation |12) 

±n (20) 
± 

where the ± denotes the boundary conditions on the two directions on the torus. The 
full partition function, given a finite discrete group G, is 

1 n (21) 

g,heG h 

where \G\ denotes the number of elements in the group. The generalization to higher 
genus Riemann surfaces is straightforward (but it is not so easy to obtain results, see 

e.g. nana). 



Z r/G -]g\ 2^ 9 



Now we specialize Eq. ([21]) to the case of the Z 2 symmetry. Since the action ( 16 1 
is invariant under g : (ft — > — 0, we have the torus partition function for the free boson 
on the orbifold gUiHIIl] 

1 



□ ++□+"□■ ( 22 ) 



y 

m " 2 

+ 

Standard CFT calculations lead to the result [T2] 

ZoM = 7i Circle ?7 + + + , 23 

2\ VdVd VdVd VdVd J 

where all the r arguments in 9 U and r\r) are understood. At the special point r\ = 1/2 
(or 77 = 2) we get 

ZorbiV = 1/2) = X 12 1 = 1 = 1 — = ^Ising • (24) 

2\ 2 I»?£>I VdVd VdVd VdVd J 8 
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Figure 2. F2(x) for the Ashkin- Teller model on the self-dual line for some values 
of rj. Inset: F%(x) — 1 in log-log scale to highlight the small x behavior. The 
black-dashed line is ~ a: 1 / 4 . 



Thus, from the orbifold partition function, using the last identity and normalizing 
such that i 7 2 4T (0) = 1, we can write the funcion F^ T (x) as 



^ 2 AT (z) = o ( F2(x) - F 2 xx (x) ) + (^(x)) 2 , (25) 



where F 2 {x) is given in Eq. (|5j), F 2 (x) is the same at rj = 1/2 (cf. Eq. (18)) 



and F 2 ls (x) is the result for Ising (cf. Eq. (111). As a consequence of the 77 O 1/rj 
symmetry of F 2 (x), also F 2 AT (x) displays the same invariance. For small x, recalling 
that F 2 (x) - 1 - tfrinlv*- 1 ] t p xx - 1 - x 1 ' 2 and F\ s - 1 ~ a; 1 / 4 , we have 




. .. for r? > 1/4, 

^'-'^ « 

The critical Ashkin- Teller model lies in the interval \/2/3 < r or b < \[2 and so 
4/3 < r\ = 2r 2 rb < 4. Thus we have F 2 AT {x) — 1 ~ a; 1 / 4 along the whole self-dual line. 
F£- T (x) for various values of r\ in the allowed range is reported in Fig. [5J where the 
behavior for small x is highlighted in the inset to show the constant 1/4 exponent. 

3. The classical Ashkin- Teller model and the Monte Carlo simulation 

The two dimensional Ashkin- Teller (AT) model on a square lattice is defined by the 
Hamiltonian 

H = 0~.[ 0~j + j'^TiTj + ViCjnTj , (27) 

(ij) (ij) 



Entanglement of 2 disjoint intervals in c — 1 theories 



where cr,; and r$ are classical Ising variables (i.e. can assume only the values ±1). Also 
the product err can be considered as an Ising variable. The model has a rich phase 
diagram whose features are reported in full details in Baxter's book [55]. We review in 
the following only the main features of this phase diagram. Under any permutation of 
the variables cr, r, err the AT model is mapped onto itself. At the level of the coupling 
constants, this implies that the model is invariant under any permutation of J, J',K. 
For K = 0, the AT model corresponds to two decoupled Ising models in a and r 
variables. For K — > oo it reduces to a single Ising model with coupling constant 
J + J' . For J = J' = K it corresponds to the four-state Potts model. It is useful to 
restrict to the symmetric Ashkin- Teller model where J = J' 

H = J) (o-iCTj + TiTj ) + K 2J 0-iCTjTiTj . (28) 

The full phase diagram is reported in Fig. [3] (in units of the inverse temperature 
/9 = 1). The model corresponds to two decoupled critical Ising models at K = and 
2J = log(l + v2)- For J = it is equivalent to a critical Ising model in the variable 
err with critical points at 2K± — ±log(l + v2). For K — > oo there are two critical 
Ising points at 2 J = ±log(l + v2). On the diagonal J — K the system corresponds 
to a 4-state Potts model which is critical at K = (log3)/4. The different kinds of 
orders appearing in the phase diagram are explained in the caption of Fig. [3j All 
the continuous lines in Fig. [3] are critical lines. The blue lines C-Is are in the Ising 
universality class. The line starting from AFIs belongs to the antiferromagnetic Ising 
universality class. On the red line ABC the system is critical and the critical exponents 
vary continuously [461 145) . 

The AT model on a planar graph can be mapped to another AT model on the 
dual graph. When specialized to the square lattice, the phase diagram is equivalent 
to its dual on the self-dual line: 

e- 2K = sinh(2 J) . (29) 

On this line, the symmetric AT model maps onto an homogeneous six-vertex model 
which is exactly solvable [45] . It follows that on the self-dual line the model is critical 
for K < (log3)/4 and its critical behavior is described by a CFT with c = 1. Along 
the self-dual line the critical exponents vary countinuously and are exactly known. 
For later convenience it is useful to parametrize the self dual line by a new parameter 
A 



4J= V252A + I e ** = l-2A, (30) 

V2 - 2A - 1 

with — 1 < A < 1/2. In terms of A, the orbifold compactification radius is [42] 
?? = 2ro2rb = 4arcc OS (-A) = 2 

where Kl is the equivalent of the Luttinger liquid parameter for the AT model. 

3.1. Cluster representation and Monte Carlo simulation 

A Swendsen-Wang type cluster algorithm for the AT model has been proposed in Ref. 
[47] and then re-derived in a simpler way by Salas and Sokal [48] . Here we partly follow 
the derivation of Salas and Sokal and we restrict to the symmetric AT Hamiltonian 



(28 1 and assume J > \K\. Using the identities for Ising type variables 

a.cij = 25 <J% „. — 1 , TiTj = 28 Ti Tj —1, (32) 
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tanh(K) 

Figure 3. Phase dia gram of the 2D symmetric Ashkin-Teller model defined by 
the Hamiltonian l ]2.s[ i . The red ABC line is the self dual line. The point B at 
K = corresponds to two uncoupled Ising models. The point C is the critical 
four-state Potts model at K = J = (log 3) /4. At J = there are two critical 
Ising points at K = ±(log(l + \f2))/2, one (Is) ferromagnetic and the other 
(AFIs) antiferromagnetic. For K — ¥ oo there is another critical Ising point at 
J = (log(l + V2))/2. All continuous lines are critical. The blue lines C — Is 
and the one starting at AFIs are in the Ising universality class. The red line 
is critical with continuously varying critical exponents. The region denoted by I 
corresponds to a ferromagnetic phase for all the variables. In the region II, cr, t, 
and err are paramagnetic. In the region III only or is ferromagnetic and in region 
IV err exhibits antiferromagnetic order while a and r are paramagnetic. 



we can rewrite Eq. ( 28 1 as 

-H= J^(2^ +2<5 TiT . -2)+Kj2( 2 ^ i a j -l){25 rtTi - 1) . (33) 

(«> (*i> 
For convenience we shift the interaction (28 1 by —4 J. In order to write the Boltzmann 
weight associated to a specific configuration we use exp(w<5 (Ti (Tj ) = (exp(w) — l)5 (Ti (T) +1 
and the analogous identity for the r variables. The Boltzmann weight of a given link 
(ij) is then 

W {ij) (vua^ruTj) = e~ 4J + [ e - 2 ( J +^) - e - 4J ] [5 ai aj + 5 n Tj ] + 

+ [1 - 2e- 2 ( J +*) + e- 4J ]^ aj S Tt Tj . (34) 

The key idea for the Swendsen-Wang algorithm is to introduce two new auxiliary 
Ising-type variables rriij and nij living on the link {ij). We redefine the Boltzmann 
weight on the link (ij) as [48] 



-4J 
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[e 2{ - J+K ^ - e AJ ][S ai<Tj 8 mij iS nij0 + S TiTj 5 mij0 6 nij i] + 

[1 - 2e- 2{J+K ^ + e- 4J ]5 CTiff A*T;<W<W ■ ( 35 ) 



Summing over my and ray we obtain the weight in Eq. (34). Eq. (J35J) has a graphical 
interpretation in terms of clusters. In fact we can divide the links of the lattice in 
"activated" (if my = 1) or "inactive" (if my = 0). The same considerations hold 
for the ray variables. Therefore, each link of the lattice can be activated by setting 
my = 1 or ny = 1. The active links connect different lattice sites forming clusters. 
There are clusters referring to the a variables (called cr-clusters) and to the r variables 
(r-clusters). Isolated lattice sites are clusters as well. Obviously, the lattice sites 
belonging to the <7-clusters (r-clusters) have the same value of a (r). The partition 



function of the extended model defined by the weight ( 35 ) can be written as 

(7,t=±1 m,n=±l (y) 

We now proceed to the following definitions. We divide all the links into three 
classes: we define Iq the total number of inactivated links; l\ the total number of 
links connecting sites which belong only to one type of clusters either a c-clustcr or a 
r-cluster. We define I2 the total number of links on which m and n are both equal to 
1. Furthermore we introduce the quantities 

B Q ^e- 4J , (37) 
Bi = [e- 2 ( J +V _ e -iJ] , (38) 
B 2 = [l- 2e- 2 ( J+K *> + e" 4J ] . (39) 



The following step is to perform the summation over a, r in Eq. ( 35 1 . This is readily 
done, obtaining the final expression for the partition function 

Z= B o ' B i B 2 2 °' T+C \ (40) 

C{r,a} 

where we denoted with C G the number of cr-clusters and with C T the total number of 
r-clusters. In the counting of r-clusters (cr-clusters) we included all the lattice sites 
connected by a link on which my = 1 (ny = 1). Isolated sites (with respect to m or 
n or both) count as single clusters. The links where my = l,ny = 1 contribute to 
both types of clusters. 



3.2. Swendsen-Wang algorithm (the direct and embedded algorithms) 

We are now in position to write the Swendsen-Wang algorithm for the symmetric AT 
model. The Monte-Carlo procedure can be divided in two steps. In the first one, 
given a configuration for (a, r) variables, we construct a configuration of the (m, n) 
variables. In the second step we update the (cr, r) variables at given (m,n). The 
details of the step one are 

• if o~i = o~j and Ti = tj, we choose (my,ny) with the following probabilities: 

- (my , ny) = (1, 1) with p 1 = 1 - 2e" 2( - J+K '> + e" 4J , 

- (my, ray) = (1,0) with p 2 = e - 2 (- 7 +^) +e~ 4J , 

- (my, ray) = (0, 1) with p 2 = e - 2( - J+K ^ + e- 4J , 

- (raiy.ray) = (0,0) with p 3 = 1- Pl - 2p 2 , 

• if <7j = o~j and r^ = —Tj, the probabilities are 

- (my, ray) = (1,0) withpi = 1 - e - 2( - J ~ K \ 
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Figure 4. A typical cluster configuration on a 12 X 12 lattice. Green lines are 
cr-clusters and red dashed lines are r-clusters. Links in blue are double links. 
Periodic boundary conditions on both directions are used. 



- (my,rijj) = (0,0) with p 2 = 1 -p u 

• if o~i = —o~j and r, = tj, the probabilities are 

- (my.ntf) = (1,0) withpi - 1 -e~ 2 ( J - K ), 

- (m ih n v] ) = (0,0) withp 2 = 1 — pi, 

• if <7j = — (Tj and r, = — Tj we choose (m.^, ny) = (0,0) with probability 1. 

In the step two, given the configuration of (to, n) generated using the rules above we 
build the connected cr-clusters and r-clusters. The value of a (r) spins are required to 
be equal within each cr-cluster (r-cluster) . We choose randomly the spin value in each 
cluster and independently of the value assumed on the other clusters. This completes 
the update scheme. (Note a typo in Ref. [25]: the minus sign in step 2 and 3 of the 
update is missing.) 

In Ref. [48] also the so called embedded version of the cluster algorithm is 
introduced. Its implementation is slightly easier compared to the direct algorithm. 
In the embedded algorithm instead of treating both a and t at the same time, one 
deals with only one variable per time. Let us consider the Boltzmann weight of a link 
(ij) at fixed configuration of r 

W {ij) (<Ti,<Tj,Ti, t 3 ) = e -^ J+K ^ + (1 - e-^ J+K ^)5 ai aj . (41) 

The model defined by this weight can be simulated with a standard Swendsen-Wang 
algorithm for the Ising model using the effective coupling constant 

■K: f -/-A--. (42) 

This is no longer translation invariant, but this does not affect the effectiveness of the 
cluster algorithm for the Ising model as long as J^j > 0. The same reasoning applies 
to the case of fixed a. Thus, the embedded algorithm is made of two steps 

• For a given configuration of r variables, we apply a standard Swendsen-Wang 
algorithm to a spins. The probability arising in the update step is pij — 

1 _ e ~2(J+KTiT 3 )_ 
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• For a given configuration of a variables, we update r with the same algorithm 
and probability = 1 - e~ 2 ^ J+KcT i a *\ 

Direct and embedded algorithms are both extremely effective procedures to sample 



the AT configurations. However, very important for the following, Eq. (40) for the 
partition function does not hold anymore for a n-sheeted Ricmann surface and we do 
not know whether it is possible to write the embedded algorithm for this case. 

3.3. Renyi entanglement entropies via Monte Carlo simulation of a classical system. 

In this section we summarize the method introduced by Caraglio and Gliozzi |18j to 
obtain the Renyi entropies via simulations of classical systems and we generalize it 
to the AT model. The partition function Z = Tr e~^ H of a (i-dimensional quantum 
system at inverse temperature /3 can be written as an Euclidean path integral in d + 1 
dimensions [5] . Thus for the n-th power of the partition function one has 

r n -f>(<M 
Z n = U V[4> k ]e *=i (43) 
J fc=i 

where <p k = <f>k(&, t) is a field living on the fc-th replica of the system and S(4>k) is 
the euclidean action (r is the imaginary time.) The actual form of the action is not 
important, but for the sake of simplicity we restrict to the case of nearest-neighbor 
interactions 

S(<t> k ) = Y,F(cp k (i), ( l> k (j)), (44) 

(ij) 

and the function F is arbitrary. We recall that Tr can be obtained by considering 
the euclidean partition function over a n-sheeted Riemann surface with branch cuts 
along the subsystem A [5]. (This equivalence is also the basis of all quantum Monte 
Carlo methods to simulate the block entanglement in any dimension |49j . ) Caraglio 
and Gliozzi constructed this n-sheeted Riemann surface for the lattice model in the 
following way. Let us consider a square lattice (for simplicity) and take the two points 
of its dual lattice surrounding A (that in 1+1 dimension is just an interval with two 
end-points). The straight line joining them defines the cut that we call A. The length 
of A is equal to the length of A. Let us consider n independent copies of this lattice 
with a cut. The n-sheeted Riemann lattice is defined by assuming that all the links of 
the fc-th replica intersecting the cut connect with the next replica fc + l(modn). To get 
the partition function over the n sheeted Riemann surface we define the corresponding 
coupled action 

n 

s n (<t> k ) = Y, E H4>k{i),4>k{3))+ F (<M*)>fc + i ( mod„)(j)). (45) 

fe=i(ij>£A (ij)ex 
This definition can be used in any dimension, even though we will use here only d = 2. 



Finally calling Z n (A) the partition function over the action (451, Tr p\ is given by 



Kpl=^- (46) 
Following Ref. [18] we introduce the observable 

Q = e -<S ra (0i,02,"-,0n;A)+2J™ =:i S(4> k :\) 
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where S n and S are the euclidean actions of the model defined on the n-sheeted lattice 
and on the n independent lattices respectively. The sum is restricted to links crossing 
the cut, as the presence of A in the arguments stresses. It then follows 



(48) 



where (•)„ stands for the average taken onto the uncoupled action Y^k=x ^i^k)- 

We can now discuss our improvement to the procedure highlighted so far. 
The practical implementation of Eq. (47 1 to calculate Trp^ is plagued by severe 



limitations: analyzing the Monte-Carlo evolution of the observable, one notices that 
it shows a huge variance because it is defined by an exponential. Direct application 
of Eq. (47) is possible then only for small lengths of the subsystem A. In order to 



overcome this problem, let us consider the quantity Z n (A) / Z n and imagine to divide 
the subsystem in L parts to have A = A\ \AA% . . . UAl, with the lengths of the various 
parts being arbitrary. Moreover we define a set of subsystems Ai = U\. =1 Ai. Then it 
holds 



Z n (A) 
Z n 



Z n {A i+ i) 



n 



(49) 



Eq. ( 49 ) is very useful because each term in the product can be simulated effectively 



using a modified version of ( 47 ) if we choose the length of Ai to be small enough. In 
fact, by definition, we have 



(°(^))^(A, 



Z n (A i+ i) 
Z n {Ai) 



where 0{Ai) is the modified observable 
We stress that in Eq 



■S n (A\)). 



(50) 



(51) 



(50 1 the expectation value in the l.h.s must be taken on t he 
coupled action on the Riemann surface with cut Ai . The disadvantage of Eq. ( 49 1 is 
that, to simulate large subsystems, one has to perform L independent simulations and 
then build the observable taking the product of the results. If the dimension of each 
piece Ai is small this task requires a large computational effort. Another important 



aspect is the estimation of the Monte Carlo error: if each term in ( 49 1 is obtained 
independently, the error in the product is 



o 



\ 



z=0 



(52) 



0{A t 



If the lengths of the intervals Ai are all equal, then the single terms of the summation 
in Eq. ( 52 ) do not change much and the total error should scale as \f~L. 

\ used another strategy to circumvent the problem with 
The trick was to consider the Fortuin-Kastelayn cluster 

The analogous for the AT 



Caraglio and Gliozzi 
the observable in Eq 



(47) 



expansion of the partition function of the Ising model 
model was reported in the previous section 

?2i oC"+C T 



Z = 



E 

C{v,t} 



D^O 7}^1 tjIi 

B Q "X ^2 1 



(53) 



where C a,T are the cr/r-cluster configurations. Going from n independent sheets to 
the n-sheeted lattice, the type of links and their total number do not change, but the 
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Tr p 0.4 




Figure 5. Trp^ for a single interval of length t in a finite system of length 
L = 120. Data have been obtained by Monte Carlo simulations using the 
embedded algorithm. The orange points correspond to the SUSY model and the 
green ones to the Z4 parafcrmions. The black crosses at I = 10 are data obtained 
using the direct algorithm. Inset: behavior of the statistical error of Trp^ vs I 
for the SUSY model. The blue-dashed line is the expected form A + Bl 1 / 2 . 



number of clusters does change, and so we get the cluster expression of observable 



(|47| for the AT model 

(D{A t ) = 2[ c -( i <+i)+ c -( i *+i)- c -( i *)-c T (^)] ( (54) 

where C (J {A i ) (C T {A i )) denote the total number of er-clustcrs (r-clusters) on the 
Riemann surface with cut A4. Since the clusters are non local objects, they represent 
"improved" observables and the variance for the Monte Carlo history of Eq. ( p34| ) is 
much smaller than in the naive implementation. 

4. The entanglement entropy in the Ashkin- Teller model 

The single interval 

We first present the results for the Ashkin- Teller model for a single interval. Although 
these results do not provide any new information about the model, they are 
fundamental checks for the effectiveness of the Monte Carlo algorithms. We performed 
simulations using both algorithms described in the previous section: the direct cluster 
algorithm and the embedded one. When using the direct algorithm, measures are 
performed using the observable |54|, while for the embedded algorithm we used the 
observable in Eq. ( 47 1 . In Fig. p3 we report the results of the simulations of Tr p\ for 
the SUSY model (r or b = \/3/2 in Fig. [l]) and for the Z4 parafermions (r or b = ^/3/2) 
both for L — 120. The orange and green points are obtained using the embedded 
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Figure 6. Plot of C2(L C ) as function of L c for different I and L. Three points 
on the self-dual line are reported: four-states Potts model, uncoupled Ising, and 



SUSY. The dashed lines are fits to the function c 2 + BL C Kl (c 2 + AL C Kl + 
BL C 2Kl for the 4-states Potts model) where is 1/2, 1, 4/3 respectively for the 



four-states Potts model, Ising, and SUSY. In the inset we report c n for n 

r 2K L /n 



for the SUSY point, 
fixed. 



The dashed lines arc fit to A + BL r 



with K l 



■ 3,4 
4/3 



algorithm. To check the implementation of the cluster observable, we report at £ = 10 
the data obtained using the direct algorithm and Eq. (54 1. The perfect agreement 
between the two results confirms the correctness of both implementations. Note that 
Tr p\ is a monotonous function of £, in contrast with the parity effects found for the 
XXZ spin chain j34j [35] that also corresponds to a vertex model [45^ . In the inset 
we show the behavior of the statistical error of the observable ( 47 1 in the SUSY case 
as function of the subsystem length I. It agrees with the prediction in Eq. ( 52 ) and 
its absolute value is extremely small, smaller than the size of the points in the main 
plot in Fig. [5] Analogous results have been obtained for all the critical points on the 
self-dual line using both algorithms. 

The results for Tr p\ in a finite system are asymptotically described by the CFT 
prediction (14 1 with n = 2 and c = 1. It is then natural to compute the ratio 

(55) 

that is expected to be asymptotically a function of the chord-length L c = [~ sin(J^)]. 
This allows to extract the non-universal quantity C2 and to check the form of the 
corrections to the scaling. In Fig. [6] we report the results for C2(L C ) for the SUSY 
point, for the two uncoupled Ising models, and for the four states Potts model. It is 
evident that for large L c , C2{L C ) approaches a constant value around 0.5. This is a 
first confirmation of the CFT predictions on the self-dual line. 
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Figure 7. P 2 ( ;E ) versus the four point ratio x for the SUSY model. The red 
points are extrapolations obtained using the finite-size ansatz ( |57| . The blue- 
dashed line is the CFT prediction. Inset: F}^ t (x) vs for the four values 
of x used in the extrap olati on (x = 0.134,0.25,0.5,0.587). The dashed lines are 
fits to finite-size ansatz ||57l. 



The previous results also provide a test for the theory of the corrections to the 
scaling to S^f . It has been shown [M] [35] that for gapless models described by 
a Luttinger liquid theory, the corrections to the scaling have the form £- 2/ W™ (or 
L c 2Kl / u for finite systems) where Kl is the Luttinger parameter, related to the circle 
compactification radius Kl — l/2rj. On the basis of general CFT arguments [3B], it 
has been argued that this scenario is valid for any CFT and so also for the AT model 
with Kl replaced by the dimension of a proper operator. It is then natural to expect 
that for the AT model this dimension is Kl in Eq. ( [31] ) , also on the basis of the results 
for the Ising model [33J [5U]. The dashed lines in Fig. [6] are fits of c 2 (L c ) with the 
function c 2 + AL~ Kl . The agreement is always very good, except for the four-state 
Potts model, for which the exponent of the leading correction Kl assumes the smallest 
value and so subleading corrections enter (as elsewhere in similar circumstances, see 
e.g. [35). In fact, the fit with the function c 2 + AL~ Kl + BL~ 2Kl is in perfect 
agreement with the data (but the presence of another fit parameter makes this result 
not so robust). This analysis confirms that Kl is the right exponent governing the 
corrections to the scaling. 

In the inset of Fig. [6] we also report the values of c n for n = 3,4 as a function of 
L c . c n becomes smaller as n increases as for the XX Z [34], XX [51], and Ising [50]l52] 
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Figure 8. F2 at (l/2) as function of r; -1 for different models (Ising, SUSY, Z4 
parafcrmions, and four-states Potts model). The blue-dashed line is the CFT 
prediction. The (colored) points close to the curve are extrapolations obtained 
with the finite-size scaling ansatz \57\ . The black crosses are the Monte Carlo 
data used for the fits. The block lengths used range from i = 5 to i = 80. 



spin-chains. The dashed lines are fits to the expected scaling behavior L c 2KL ' n of the 
corrections, that reproduce perfectly the data. 

4-2. The entanglement entropy of two disjoint intervals. 

In this section we investigate the entanglement entropy of two disjoint intervals and 



check the correctness of our prediction (25) for the AT model on the self-dual line. 



As for all other cases studied so far numerically (i.e. Heisenberg [17], Ising [24l 126] . 
and XY [26] chains) , strong scaling corrections affect the determination of the scaling 
function F n (x). CFT predictions have been confirmed only using the general theory 
of corrections to the scaling [331 [SHI |3_7] . 

In order to determine the function F n (x), we consider the ratio 

Tr o n 

F - t{x) = Trn» Tr^» (1 " , (56) 

lr PAi 1t Pa 2 

and, on the basis of the general CFT arguments [3B], we expect that the the leading 
correction to scaling can be effectively taken into account by the scaling ansatz 

F^(x) = FZ FT (x) + r 2 ^ n f n (x) + .... (57) 

For the Ising model it has been found to = 1/2 [231 12S]- Since for 77 = 2 the AT 
Hamiltonian reduces to two uncoupled Ising models, one naively expects u) — 
along the whole self-dual critical line of the AT model. 

Hereafter we only consider Tr p\- We start our analysis from the SUSY point 
that (assuming ui = should have the smaller corrections to scaling. In Fig. [7] 

we show Monte Carlo data at I = 10,20 (L = 120) for F^x) plotted against the 



four point ratio x defined as in Eq. (15). We report with the blue dashed line the 
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Figure 9. F 2 CFT (l/2) - Fj at (l/2) versus The dashed lines are fits to the 

finite-size scaling ansatz 
in log-log scale. 



|57| fixing the value of F 2 CFT (l/2). Left: the same plot 



asymptotic CFT result (cf. Eq. ( 25 )). As in all other cases considered in the literature 
[Ml the curves for Fi(x) at I = 10, 20 are not symmetric functions of x — >• 1 — x, 
as instead the asymptotic CFT prediction must always be [IT]. This is due to the 
non-symmetrical finite-size corrections fa(x) in Eq. (57 1. We extrapolate the result 
at £ — > oo using the ansatz (57) and u> = 2/3. The extrapolations are reported as red 
points in Fig. [7] There is a very good agreement between the extrapolations and the 
theoretical curve. Since the correction exponent lu = 2/3 is rather large, and so the 
corrections small, even small subsystems such as I — 10, 20 are enough to obtain a 
good extrapolation. In the inset of Fig. [7] we report the Monte Carlo data for Fl &t (x) 
against t~ 2 l z . The linear behavior in this inset confirms the validity of the ansatz (57) 
and the reported straight lines are the fits giving the extrapolations reported in the 
main panel. 

We also investigate other points on the self dual line, namely the 4-states Potts 
model (77 = 4), the parafermion Z4 (n = 3), the uncoupled Isings (j] = 2). In Fig. [8] we 
report F^(x) - Ff? FT (x) at fixed x = 1/2 versus rj for all the mentioned models, 
we report x = 1/2 because it is the value of x providing the most stable estimate, 
but also other values have been studied. Indeed, on one hand, the computational cost 
of the simulations decreases going toward x — 1 (the reason being evident from the 
definition of x for which smaller lattice sizes are needed). On the other hand, scaling 
corrections become more severe in the region x ~ 1, as clear from the results for the 
SUSY model in Fig. [7] Thus x = 1/2 represents the best compromise between these 
two drawbacks. The dashed curves in the left panel of Fig. [9] are fits of the data 
with Eq. (57) obtained by fixing the value of F^^ix) to its predicted value (cf. Eq. 
(25)). There is a very good agreement with the full theoretical picture, confirming 
in particular the correctness of the exponent governing the leading correction to the 
scaling. For the Z4 parafermions and for the four-state Potts model, we needed very 
large values of £ in order to show the correct asymptotic behavior (the range of £ 
reported in the plot is in fact 5 < £ < 80). This is made clearer in the right panel 
of Fig. [9] where the same data are shown in log-log scale. In Fig. [8] we reports the 
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Figure 10. Monte Carlo data for f2(x) obtained as $i(x) = (F^ix) — 
F§ FT (x))£ Kl ^ 2 as function of x. We show data for £ = 10 and various 
models (SUSY, Z4 parafermions, Ising model and the model corresponding to 
rf~ 1 = 0.74). The blue-dashed lines are asymptotic fits to Ax 1 / 4 . 



fits obtained by fixing only the exponent of the corrections oj = Kl/2 and leaving 
F^ FT (l/2) free. For all considered values of 77, the extrapolation of F2 at (l/2) to 
t — >• 00 is compatible (within error bars) with the expected result F^ FT (l/2). 

We finally study the correction amplitude fa{x) in Eq. ( |57[ ). This function is the 
main reason of the asymmetry in x — > 1 — x for F2 (%) and knowing its gross features 
could greatly simplify future analyses. For the Ising model, it has been found that 
fz(x) ~ x 1 / 4 for small x, that is the same behavior of F<z(x) — 1. Since along the whole 
self-dual line ^(ar) — 1 ~ 2 1 / 4 , we would expect 

/a(x)~a! 1/4 . (58) 

For the Ising model (i.e. 77 = 2), this scenario has been already verified with high 
precision |24|_ 



In Fig. 10 we report fa{x) obtained as f 2 (x) = (F 2 CFT (x) -F^ t )i KL l 2 as function 
of x (in logarithmic scale to highlight the small x behavior) . All data correspond to 
£ = 10 and various values of L. For the two largest values of 77 (Z4 parafermionic 
theory at 77 = 3 and for the Ising model at 77 = 2), we observe an excellent agreement 
with our conjecture f2(x) ~ x 1 ^ 4 . However decreasing the value of 77, i.e. for the 
SUSY model at 77 = 3/2 and for the model at n^ 1 — 0.74, the behavior of f2(x) is not 
as linear as before, especially for high value of x. Nonetheless for x < 0.4 the data 
confirm the behavior x 1 ^ 4 . Furthermore, it seems that for any 77 ^ 2, subleading terms 
in the expansion for small x appear and they are vanishing only for the Ising model. 

5. The Tree Tensor Network 



This section is divided into two parts. First we explain in a self contained way how 
to extract the spectrum of the reduced density matrix of some specific bipartitions 




of a pure state encoded in a Tree Tensor Network (TTN). We only recall the 
basic definitions introduced in Ref. [S3j and refer the reader to the literature for 
complementary works on the subject [SlEiEiETlEllSaiSnilSIlEaiMlIMlEaiM]^ 
Secondly we quickly recall how to use TTN to calculate the ground state of the 
anisotropic Heisenberg spin-chain. 



5.1. Tree Tensor network and reduced density matrices. 

We consider a one dimensional lattice C made of N sites, where each site is described 
by a local Hilbert space V of finite dimension d. In this work the state is the ground 
state I^gs) of some local Hamiltonian H defined on £, but in general it could be an 
arbitrary pure state \^>) E Y® N defined on the lattice C. 
A generic state € ~V® N can always be expanded as 

d d d 

W=E E ••• E Ti li2 ... iff \k)\i 2 )---\i N ), (59) 

tl = l 12 = 1 i]V = l 

where the d N coefficients Ti 1 i 2 ...i N are complex numbers and the vectors 
|2 S ), • • • , \d s )} denote a local basis on the site s E C. We refer to the index 
i s that labels a local basis for site s (i s = 1, • • • , d) as a physical index. 



In the case we are interested in, the tensor of coefficients T ili2 ... iN in Eq. (59) 
is the result of the contraction of a TTN. As shown in Fig. [TT] for lattices of = 4 
and N = 8 sites, a TTN decomposition of T ili2 ... iN consists of a collection of tensors 
w that have both bond indices and physical indices. The tensors are interconnected 
by the bond indices according to a tree pattern. The N physical indices correspond 
to the leaves of the tree. Upon summing over all the bond indices, the TTN produces 
the d N complex coefficients T ili2 ... iN of Eq. (59). 
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Figure 12. (i) Diagrammatic represent atio n of the two types of isometric tensors 
in the TTN for a N = 4 lattice in Fig. |11| (ii) Graphical representation of the 
constraints in Eqs. (161} and K2h fulfilled by the isometric tensors. 



The tensors in the TTN will be constrained to be isometric, in the following sense. 
As shown in Fig. [12] for the N = 4 lattice of Fig. [IT] each tensor w in a TTN has 
at most one upper leg/index a and two lower indices/legs Px,^2, so that its entries 
read (w)^ ^ 2 (everything can be generalized to tensors with more upper and lower 
legs [53]). Then we impose that 

EH^> f )*— ( 6 °) 

ft ,ft 

For clarity, throughout this paper we use diagrams to represent tensors networks as 
well as tensor manipulations. For instance, the constraints for the tensors w\ and W2 
of the TTN of Fig. [TTJfor a N = 4 lattice, namely 

EMftft^)"—^'' ( 61 ) 

ft ft 

=1- (62) 

ft ft 

are represented as the diagrams in Fig. [l2|ii) . We refer to a tensor w that fulfills Eq. 
( 60 1 as an isometry. 

An intuitive interpretation of the use of a TTN to represent a state I*) can be 
obtained in terms of a coarse-graining transformation for the lattice C. Notice that 



the isometries w in Fig. 11 are organized in layers. The bond indices between two 



layers can be interpreted as defining the sites of an effective lattice. In other words, 
the TTN defines a sequence of increasingly coarser lattices {Cq, Ci, ■ ■ ■ , Ct-i}, where 
£ = £ and each site of lattice C T is defined in terms of several sites of £ r _i by means 
of an isometry w T , see Fig. [13] In this picture, a site of the lattice C T effectively 
corresponds to some number n T of sites of the original lattice Cq. For instance, each 



of the two sites of Li in Fig. 13 corresponds to 8 sites of Co- Similarly, each site of 
lattice Ci corresponds to 4 sites of £q. 

The use of isometric tensors, and the fact that each bond unambiguously defines 
two parts (A : B) of the chain which are connected only through that bond as displayed 
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Figure 13. The isometric TTN of Fig. [TTJfor a N = 8 lattice Co with periodic 
boundary conditions (the blue external circle) is associated with a coarse-graining 
transformation that generates a sequence of increasingly coarse-grained lattices 
C\, C2 and £3 (the inner circles). Notice that in this example we have added an 
extra index to the top isometry 103, corresponding to the single site of an extra 
top lattice £3, which we can use to encode in the TTN a whole subspace of ~V® N 
instead of a single state 



in Fig. [TJJ implies that the rank of that bond in the TTN is given by the Schmidt 
rank \(A : B) of the partition (A : B) [53]. Thus the reduced density matrix pa for 
a set A of sites of C is 

p A =tr B |*X*l=EP-l*«X*«l' ( 63 ) 

a 

where p a are the eigenvalues of pa- It follows then the Renyi entanglement entropies 

Cl( n ) 

S A are 

a 

and for n = 1 

Sa ] = -tr(p A logp A ) = -J^PalogPa- (65) 

a 

In the following we denote the ranks of the tensor w T , a, Pi, /3 2 as , \ T , x T ~ ,X T ~ ■ 
In general, they fulfill 

X T < (X^ 1 ) 2 , (66) 
meaning that w T projects states in V r_1 ® V r_1 into the smaller Hilbert space V T . 




Figure 14. By erasing one of the indices in the TTN the spin chain is always 
divided in two parts A and B |59l . Here we show that in the case of the N = 8 
lattice of Fig. |ll| there are three classes of indices, identified by their position in 
the TTN. i) physical bonds connect a single spin with the rest of the lattice, ii) 
bond indices of the first layer connect a block of two adjacent spins to the rest of 
the lattice, iii) bond indices of the third layer of the lattice connect four adjacent 
spins, to the other half. This implies that the rank of the index is the Schmidt 
rank of the respective partition. 



For a critical chain, the logarithmic scaling of the entanglement entropy (cf. Eq. 
([2])) implies that the rank of the isometries should at least grow proportionally to the 
length of the block represented by the effective spins 

X T oc n T , (67) 

which means that while moving to higher layer of the tensor network the rank of 
the isometries increases. This also implies that the leading cost of the computation is 
concentrated in contracting the first few layers of the TTN. li N = 2 T and we describe 
a pure state (so that the rank of the a T is one) the maximal rank of the tensors in the 
TTN is 

X = maxx r = X T_1 . (68) 

T 

In Ref. [53J it has been shown that i) a TTN description of the ground state of 
chain of length N with periodic boundary conditions can be obtained numerically with 
a cost of order C(log Nx 4 )- ii) From the TTN it is also straightforward to compute 
the spectrum {p a } of the reduced density matrix pA (cf. Eq. (63)) when A is a 
block of contiguous sites corresponding to an effective site of any of the coarse-grained 



lattices Ci, ■ • • , £t-i- Fig. 15 illustrates the tensor network corresponding to pa for 



the case when A is one half of the chain. Many pairs of isometries are annihilated. 
In addition, the isometries contained within region A can be removed since they do 
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Figure 15. Computation of the spectrum {p a } of the reduced density matrix 
pA for a block A that corresponds to one of the coarse-grained sites, (i) Tensor 
network corresponding to pa where A is half of the lattice, (ii) Tensor network 
left after several isometries are annihilated with their Hermitian conjugate, (iii) 
since the spectrum of pa is not changed by the isometries acting on A, we can 
eliminate them and we are left with a network consisting of only two tensors, which 
can now be contracted together. The cost of this computation is proportional to 

0(x'V) < c(x 5 )- 



not affect the spectrum of pa- From the spectrum {p a }, we can now obtain the Renyi 

(n) 

entropies S A . The leading cost for computing the spectrum of the reduced density 
matrix pa for this class of bipartitions is due to the contractions of the first layers of 
the TTN. When the bipartition is such that A is a quarter of the chain, this implies 
a cost proportional to 0(x' 3 X 2 ) < x'x 4 > where x' = X T ~ 2 - 

It is also possible to compute the reduced density matrix pa when A is composed 
of two disjoint subintervals A\ and A2, where now each of the two intervals is a block 
of contiguous sites corresponding to an effective site of the coarse grained lattice. The 
cost of this computation is again dominated by contracting the upper part of the tensor 
network, and the most expensive case is obtained by considering A as the collection of 
two N/4 spins blocks, separated by N/4 spins. The tensor network corresponding to 



this pa is shown in Fig. 16 Also in this case many pairs of isometries are annihilated. 
The isometries contained within the composed region A can also be removed since 
they do not affect the spectrum of pa- The cost of contracting this tensor network is 
proportional to max[C(x 2 x' 4 ), C(x 3 x' 2 )] < C(x 6 )- 



5.2. The TTN and the anisotropic Heisenberg spin-chain 

In the previous subsection we have shown how to extract the spectrum of the reduced 
density matrix for a single and a double spin block from a TTN state. In this 
manuscript we are interested in reduced density matrices calculated on the ground- 
state of the anisotropic Heisenberg spin chain (XXZ model) in zero magnetic field, 
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3ZH 



in 



iii) 



Figure 16. Computation of the spectrum {p a } of the reduced density matrix 
PA when A corresponds to two coarse-grained sites separated by one coarse 
grained site from both sides, (i) Tensor network corresponding to pa where 
A is a quarter of the lattice, (ii) Tensor network left after several isometries 
are annihilated with their Hermitian conjugate, (iii) Since the spectrum of 
pA is not changed by the isometries acting on A, we can eliminate them and 
we are left with a network consisting of only few tensors, which can now be 
contracted together. The cost of contracting this tensor network is proportional 
to max[0(xV 4 ),0(x 3 x' 2 )] < 0(x 6 )- 



defined by the Hamiltonian 

H = X>K+i + W+i + ■ ( 69 ) 

3=1 

where a" are the Pauli matrices at the site j. Periodic boundary conditions are 
assumed. We are interested in gapless conformal phases of the model, that is 
— 1 < A < 1. This phase is described by a free-bosonic CFT compactified on a 
circle with radius that depends on the parameter A 

9 1 arccosf— A) . , 

V = 2r c 2 irclc = wjt = " , (70) 

ZK L 7T 

where K L is the Luttingcr liquid parameter. [§] The sign convention in the Hamiltonian 



(69 1 is such that the model is (anti)ferromagnetic for A < (A > 0). Hamiltonian 
(69 1 is diagonalizable by means of Bethe ansatz. However, obtaining the spectrum of 
the reduced density matrix from Bethe ansatz is still a major problem and only results 
for small subsystems are known [44j [67] . For this reason we exploit variational TTN 
techniques to obtain the ground state. 

Here we follow the variational procedure described in detail in Ref. |53j . where 
the generic technique (consisting of assuming a tensor network description of the 
ground state and minimize the energy variationally improving the tensors one by one 
as described, i.e., in Ref. [2]) has been specialized and optimized for the case of a 

§ Notice similarities and differences between Eq. j70| l and its analogous for the AT model ( |3 1 [ ) . The 
relation between r\ and r 2 and the relation between Kl and A are the same for both XXZ spin-chain 
and AT model, but the relation between 77 and A (or and r) is different. 
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Figure 17. TTN data for the non universal constant C2(L C ) as function of the 
chord length L c for different values of A. The dashed curves are fits to the 
function A + BL C Kh . The reported data have been obtained with L = 128 for 
A = 0, 0.1, 0.6 and L = 64 for the other values. 



TTN. We exploit translation invariance by using the same tensor at each layer of the 
TTN. One could also improve the efficiency further by exploiting the U{X) symmetry 
of the Hamiltonian (69), i.e. the rotations around the z axis. However we did not 
make use of this symmetry here. 



6. The Block Entanglement of the Anisotropic Heisenberg spin-chain 

In this section we report the TTN results for the Renyi entropies in the XXZ spin- 
chain for a single and a double interval. As a main advantage compared to the classical 
Monte Carlo simulations performed for the AT model, with a single TTN simulation 
we obtain the spectrum of the reduced density matrix and hence any Renyi entropy, 
including von Neumann Oppositely with the Monte Carlo methods only Renyi 

(n) 

entropies S\ of integer order n > 2 can be obtained and each of them requires an 
independent simulation. 



6.1. The single interval. 

We first present the TTN results for the single interval. These have been already 
obtained with many numerical variational techniques [34 , 68 , 69 , 38J and are reported 
here only to test the accuracy of the TTN and to fix units/scales etc. Using variational 
TTN, we find the ground-state of the XXZ Hamiltonian ( [69] ) and from this we extract 
the spectrum of the reduced density matrix of the single block, as explained in the 
previous section. We then numerically obtain Tr p n A . The maximum size of the chain 
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Figure 18. TTN data for F^ it (x) as function of x for various sizes of the chain 
L = 16, 32, 64, 128, subsystem lengths i = 4, 8, 16, 32, and A = -0.3, -0.1, 01, 0.6. 
Different values of A are distinguished by different colors, while different symbols 
denote different values of i. The arrows denote the (asymptotically) increasing 
subsystem sizes I. 



that we consider is L = 128. The subsystem lengths considered are I = 2, 4, 8, 16, 32. 
Notice that with the TTN method, using a binary tree as we are doing, we can 
effectively access only subsystems sizes of the form 2 m with m arbitrary integer, as it 
should be clear from the previous section. In particular this limits the calculation to 
even values of £ and we can not study the parity effects reported in Ref. [3H [35] . 

We considered different values of the anisotropy parameter A, namely A = 
-0.3,-0.1,0,0.1,0.2,0.4,0.6,0.8,1. The TTN becomes less effective for values of 
A < —0.5. This can be easily traced back to the smallness of the finite-size gap that in 
the minimization process causes the algorithm to be stuck in meta-stable states when 
the system size is large enough. This drawback could be cured by using larger values 
of x (and so larger computational cost), but as we shall see, the considered values 
of A suffice to draw a very general picture of the entanglement. For the isotropic 
Heisenberg antiferromagnet at A = 1 we ignore the presence of logarithmic corrections 
to the scaling [551 [Mj; that have a minimal effect for all our aims. 

As for the AT model, we study the quantity C2(L C ) defined by the ratio in Eq. 



(55 1. The results are shown in Fig. 17 for all considered values of A. The scaling 



corrections arc evident, especially for larger values of A, as expected [33]. These 
corrections for Trp^ are indeed of the form L c 2Kh ^ n [jg] (K^ is defined in Eq. (70)). 
The dashed lines reported in Fig. [IT] are fits to this form for n = 2, showing the 
agreement between TTN data and the fits. We checked that all the TTN data agree 
with the ones obtained in Ref. 34J using density matrix renormalization group. The 
agreement is perfect and for this reason we refer to the above paper for a detailed 
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Figure 19. TTN data for F^ at (l/2) 



F 2 CFT (l/2) as function of l/l for various 



A. The dashed lines are fits to the function with the generalized finite- 

1ZD- 



ansatz 



study of Trp^ for n > 2. 



6.2. Double interval: the n = 2 case. 

We now consider a subsystem made of two parts A\ and of equal length I. We 
start by studying the quantity Tr p\ iU a 2 for finite chains and extract the universal 
function T (:r) by proper extrapolation. Since we only consider even £, corrections 
to the scaling are expected to be monotonic in I also for F?,{x), oppositely to the case 
of arbitrary I parity [T71 HB] ■ Tbe CFT prediction for the function F 2 (x) for the XXZ 
chain is Eq. j|5| with rj given by Eq. ( 70 ) . 

In Fig. [T8] we report TTN data for F^ (x) (obtained with the ratio defined in 
Eq. (56 1) as function of the cross ratio x for A = —0.3,-0.1,0.1,0.6 and subsystem 
sizes £ = 4,8, 16, 32. The different values of A are denoted with different colors, while 
the different symbols stand for the various I. On the same figure we also show the 
asymptotic F^ FT (a:) as dashed lines. It is evident that strong scaling corrections 
affect the data, as expected. Colored arrows denote the direction of (asymptotically) 
increasing subsystem sizes. Very surprisingly, while for A = —0.3,-0.1,0.1 the 
asymptotic CFT result is approached from below, for A = 0.6 it is approached from 
above. Moreover, for A = 0.6 the behavior of the data is not monotonic. This 
contrasts the results obtained for the AT model in the previous sections and the ones 
obtained for the XX and Ising spin-chains [26 . 

In order to shed some light on this unexpected pheno men on, it is worth to look 

we report one of these 
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at i 7 ^*^) as functions of I for fixed values of x. In Fig 
plots for x = 1/2. Analogous figures are obtained for other values of x. Corrections 
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Figure 20. TTN data for F^(x) as function of x for various sizes of the chain, 
A = -0.3,0.1,0.6,1, and subsystem lengths t = 4,8,16,32. We denote with 
different symbols the values of i and with different colors the various A. The 
dashed curves are the theoretical results given by Eq. The arrows denote the 
(asymptotically) increasing subsystem sizes I. 



to the scaling are non-monotonic in the range 0.2 < A < 0.7. This phenomenon can 
be understood if further corrections to the scaling are taken into account. There are 
two corrections that can be responsible of this behavior. On the one hand, corrections 
of the form £~ mK L (from £- 2mK ^/n at n = 2) for any integer m are know to be 
present [3S], on the other hand usual analytic corrections such as l~ x are generically 
expected to exist for any quantity from general scaling arguments. Thus the most 
general finite-^ ansatz has the form 

Ft\ X ) = F^( X) + M^i + m + i0i . . . , (71) 

where the first correction is the unusual one employed also for the Ashkin- Teller model, 
and the other two are the ones just discussed. The effect of subleading corrections 
is enhanced by the fact the the amplitude functions fi{x) and Ja{x) or /b(x) have 
opposite signs determining the non-monotonic behavior. Unfortunately, for values of A 
for which the effect of subleading corrections is more pronounced (i.e. 0.1 < A < 0.6), 
wc have K L < 1 < 2K L , making difficult to disentangle corrections with close 
exponents. Thus, in order to present analyses of a good quality, we ignore the last 
correction (i.e. we fix /s(a;) = 0). To check the proposed scenario, we performed the 
fit of the data in Fig. [l9]with the ansatz (71) and = 0. The results of the fits 

are reported in the same figure, showing perfect agreement with the data for all the 
values of A. We repeated the same analysis for other values of x, finding the same 
quality of fits as for x = 1/2. However, we cannot exclude that corrections of the form 
£-2K L nave an important role. 
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Figure 21. TTN data for F' at (:r) at fixed x = 
I = 4,8, 16,32, The considered values of A are A 
dashed curves are fits with the ansatz 1721 



1/2 as function of I 1 for 
= -0.3, -0.1,01,0.6, 1. The 



6.3. Double interval: the n = 3 case. 



Now we report the same analysis performed for Tr p\ for the third moment of pA, 
i.e. Tr p\. Again we consider finite-size XXZ spin-chains and extract the universal 
function Fpf T (x) by finite-size analysis. The expected CFT result is given for general 
n by Eq. Q. In Fig.|o]we show TTN data for F^(x) (obtained from Eq. @) at 



A = —0.3, 0.1, 0.6, 1 and subsystem sizes up to I = 32. We also show the theoretical 
curves given by Eq. Q. As for the n — 2, the asymptotic universal curve is approached 
from below for A < 0.6, and from above for A > 0.6. Furthermore, the behavior of 
the numerical data for A > 0.6 is non monotonia This suggests that the ansatz in 



Eq. (57) is not enough to describe accurately the TTN data and further corrections 
to the scaling should be included as for Trp^. 

For n — 3, the leading corrections to the scaling are described by the ansatz ([57]), 
i.e. the leading exponent is 2Kl /3. Thus, for the cases when subleading corrections are 
more important (i.e. for A > 0.6) the ordering of the exponents is 2K^/3 < < 1 

and so it is reasonable to ignore the analytic correction. Thus we fit TTN data with 
the function 

Ft\x) - F 3 CFT (x) = / 3 (x)r 2 **/ a + f B {x)r iK -'* . (72) 



In Fig. 21 we report TTN data for F^ t {x) - F 3 ^ FT (x) for 1=1/2 and several values 



of A. The dashed lines are fits with the finite-size ansatz ( |72| ), that perfectly reproduce 
the data. 

6.4- Double interval: The von Neumann entropy. 

TTN gives access to the full spectrum of the reduced density matrix of A\ U A-i and so 
to the entanglement entropy as well. In Fig. 22 we report the function Fy^ r (x) 
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Figure 22. TTN data for the von Neumann entropy for various values of A 
in the interval [—0.3, 1]. We show with different symbols the values of A while 
different colors stand for different £ and lattice sizes. 



defined as 

F%(x) = S$ uAa S$ ~ S<£ - i log(l - x) , (73) 

for A in the interval [—0.3, 1] for various L up to 128 and subsystem sizes I — 
2,4,8,16,32,64. We indicate with different symbols different values of A, while the 
colors are for various sizes t. As known from many other investigations on single 
and double intervals (quantum Ising spin chain, XY model, XXZ) the von Neumann 
entropy does not show oscillations with the parity of the subsystem and the corrections 
are much smaller, actually negligible from any practical porpouse. Fig. [22] confirms 
this observation for the two interval entanglement entropy for the XXZ spin-chain in 
a wide range of A. Indeed, at fixed value of A perfect data collapse is observed even 
for very small values of £. 

Unfortunately, as already stated in the introduction, the CFT prediction for 
Fvn{x) is unknown for general x because the analytic continuation of F n (x) to non- 
integer n is not achievable. However, an expression for the leading term of the small 
x expansion of Fvn(x) has been recently extracted [3D] from Eq. (13) 

f vn ( X ) = [- a ) ^r (Q + 3/V (74) 

where a = min[ry, I/77] and T is the Euler function (not to be confused with the T 
matrix in Eq. pi). In order to check the correctness of this formula, in Fig. 23 we 



report the same data for Fvn(x) in a log-log scale to highlight the power-law behavior 
for small x. We also report the small x expected from Eq. (74 1. For A = —0.3 the 
agreement is good, but it gets worse increasing A. The natural explanation is that 




Figure 23. TTN data for the von Neumann entropy for A = -0.3,0.1,0.6,0.8 
(the data for different A are denoted with different symbols) in log-log scale. We 
used different colors to indicate the different block sizes i and lattice sizes L. The 
continuous lines are the small x behavior obtained from \74\ . The dashed lines 
are the small x behavior where the O(x) term has been added as in Eq. {75}. 



the considered values of x are not small enough for the asymptotic Eq. ( 74 ) to be 
valid. We should then include further terms in the small x expansion. As explained 
in Ref. |30j . further coefficients in the expansion for small x are difficult to obtain in 
general. However, there is a term that is very easy to obtain and that (luckily enough) 
is responsible of the previous disagreement. Indeed, as shown in Ref. [30] (cf. Eq. 70 
and 71 there) the function F n (x) has always (i.e. independently of 77) a simple 0(x) 
contribution coming from the denominator in Eq. Q, i.e. |6(0|r)| 2 = l+x(n-l/n)/6, 
that can be easily analytically continued giving 

x \° nz r (« + !) x , mJioL\ 



Fvs{*)=\j) ^^-3+0(0- (75) 

Notice that the added term becomes more important when a is close to 1, i.e. in the 
XXZ spin-chain when A approaches 1. In Fig. 23 we also report the prediction (75) 
as dashed line, that is asymptotically in perfect agreement with the numerical data 
for all values of A. 



7. Conclusions 

In this manuscript we provided a number of results for the asymptotic scaling of the 
Renyi entanglement entropies in strongly interacting lattice models described by CFTs 
with c = 1. Schematically our results can be summarized as follows. 

• We provided the analytic CFT result for the scaling function F2(x) for S A in 
the case of a free boson compactified on an orbifold describing, among the other 
things, the scaling limit of the Ashkin- Teller model on the self-dual line. The 
final result is given in Eq. ( 25 ) . 
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We developed a cluster Monte Carlo algorithm for the two-dimensional Ashkin- 
Teller model (generalizing the procedure of Caraglio and Gliozzi [TH] for the 
Ising model) that gives the scaling functions of the Renyi entanglement entropy 
(for integer n) of the corresponding one-dimensional quantum model. With this 
algorithm, we calculated numerically the scaling function ^2(2;) of the AT model 
along the self-dual line and we confirm the validity of the CFT prediction. In 
order to obtain a quantitative agreement, the corrections to scaling induced by 
the finite length of the blocks are properly taken into account. 

We considered the XXZ spin chains by means of a tree tensor network (TTN) 
algorithm. The low-energy excitations of model are described by a free boson 
compactified on a circle for which CFT predictions are already available both for 
n = 2 [T7] and for general integer n 19J . Taking into account the corrections to 
the scaling, we confirm these predictions (that resisted until now to quantitative 
tests) for n = 2,3. Furthermore, we provide numerical determinations of the 



scaling function of the von Neumann entropy (cf. Fig. 22 1 for which CFT 
predictions do not exist yet for general x. For small x we confirm the recent 
prediction of Ref. [3(1 (cf. Fig. ||). 

The methods we employed (classical Monte Carlo with cluster observables and 
TTN) are very general techniques that can be easily adapted to other models of 
physical interest. On the CFT side, it must be mentioned that a closed form for the 
functions F n (x) at integer n for a free boson compactified on an orbifold is not yet 
available, but work in this direction is in progress [70]. 
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